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ABSTRACT 


The objective of this proposed project is to research and develop a prediction tool for advanced 
additive manufacturing (AAM) processes for advanced materials and develop experimental 
methods to provide fundamental properties and establish validation data. Aircraft structures and 
engines demand materials that are stronger, useable at much higher temperatures, provide less 
acoustic transmission, and enable more aeroelastic tailoring than those currently used. Significant 
improvements in properties can only be achieved by processing the materials under non- 
equilibrium conditions, such as AAM processes. AAM processes encompass a class of processes 
that use a focused heat source to create a melt pool on a substrate. Examples include Electron 
Beam Ereeform Eabrication and Direct Metal Deposition. These types of additive processes 
enable fabrication of parts directly from CAD drawings. 

To achieve the desired material properties and geometries of the final structure, assessing the 
impact of process parameters and predicting optimized conditions with numerical modeling as an 
effective prediction tool is necessary. The targets for the processing are multiple and at different 
spatial scales, and the physical phenomena associated occur in multiphysics and multiscale. 

In this project, the research work has been developed to model AAM processes in a multiscale 
and multiphysics approach. A macroscale model was developed to investigate the residual 
stresses and distortion in AAM processes. A sequentially coupled, thermomechanical, finite 
element model was developed and validated experimentally. The results showed the temperature 
distribution, residual stress, and deformation within the formed deposits and substrates. 

A mesoscale model was developed to include heat transfer, phase change with mushy zone, 
incompressible free surface flow, solute redistribution, and surface tension. Because of excessive 
computing time needed, a parallel computing approach was also tested. In addition, after 
investigating various methods, a Smoothed Particle Hydrodynamics Model (SPH Model) was 
developed to model wire feeding process. Its computational efficiency and simple architecture 
makes it more robust and flexible than other models. More research on material properties may be 
needed to realistically model the AAM processes. 

A microscale model was developed to investigate heterogeneous nucleation, dendritic grain 
growth, epitaxial growth of columnar grains, columnar-to-equiaxed transition, grain transport in 
melt, and other properties. The orientations of the columnar grains were almost perpendicular to 
the laser motion’s direction. Compared to the similar studies in the literature, the multiple grain 
morphology modeling result is in the same order of magnitude as optical morphologies in the 
experiment. 

Experimental work was conducted to validate different models. An infrared camera was 
incorporated as a process monitoring and validating tool to identify the solidus and mushy zones 
during deposition. The images were successfully processed to identify these regions. 

This research project has investigated multiscale and multiphysics of the complex AAM 
processes thus leading to advanced understanding of these processes. The project has also 
developed several modeling tools and experimental validation tools that will be very critical in 
the future of AAM process qualification and certification. 
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1. Introduction 


Table 1.1 lists the tasks completed for the entire project. In the main body, chapter 2 describes the 
three-dimensional model used to simulate the thermal and velocity field of the laser deposited 
melt pool. Chapter 3 introduces the idea for the wire model, which solves the solid-liquid 
interaction and shows preliminary simulation results for the wire and molten pool. Chapter 4 
describes the use of parallel computing in order to optimize the thermofluid simulation previously 
mentioned in chapter 2. Chapters 5 and 6 introduce the Smoothed Particle Hydrodynamics Model 
(SPH Model), a comparison with the Volume of Fluid (VOF) technique, and the SPH Model 
implementation. Chapter 7 describes the numerical analysis of thermal stress and deformation in 
multilayer Laser Metal Deposition Processes (LMDP) and the experimental validation. Chapter 8 
covers the Finite Element-Cellular Automaton (FE-CA) model and predicts equiaxed and 
columnar grain growth based on the thermal history obtained with the computations described in 
Chapter 7. The grain growth prediction was validated by comparing similar experimental results 
from the literature. Chapter 9 introduces the microstructural observation and simulation of 
Ti6A14V single- wall deposition. Chapter 10 is about the use of an infrared camera to monitor the 
laser metal deposition process and the concomitant image processing. 
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Table 1.1 Tasks of the project 


Task 1. macroscopic model 


Task la. Thermal/fluid dynamics 

Code completed/Python version developed to 
improve code efficiency and numerical 
convergence issues 

Task lb. Residual stress model 

Completed using Abaqus commercial software 

Task 2. Microscopic model 


Task 2a. As-received microstructure 

Completed using open source MMSP code 

Task 2b. Solidification microstructure 

2D code completed. Columnar model completed; 
3D model is in progress. 

Task 2c. Solid-state phase transformations 

Being implemented 

Task 3. Data collection and model validation 


Task 3a. Temperature measurement 

Technique developed (IR camera) and thermal 
couple to be used for validation 

Task 3b. Clad and melt pool geometry 

Technique developed; used for validation 

measurement 


Task 3c. Residual stress measurement 

Completed using a simple example to demonstrate 

Task 3d. Validation of microstructures 

LAMP experiments completed. 

Task 3e. Validation using NASA's EBF^ apparatus 

In progress; wire is being modeled in the 
code, initial result has been shown 
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2. Thermo-Fluid Model 


A 3D numerical model for laser deposition by powder injection was developed, and simulation 
codes were run for different process parameters. In this model, the Piecewise Linear Interface 
Calculation (PLIC) scheme [1][2][3], which approximates the interface as a plane within a 
computational cell, was used to track the free surface. Compared to earlier versions of the 
Volume Of Fluid (VOF) methods [4] [5] [6] [7], PLIC-VOF is more accurate and avoids the 
numerical instability [8]. It also can handle large interface deformation. Surface tension was 
modeled as a volume force using the Continuum Surface Force (CSF) model [9], combined with a 
force-balance flow algorithm [10]. The basic idea underlying the CSF model is the representation 
of surface tension as a continuous force per unit volume that acts in a finite transition region 
between the liquid and gas phases. The interface curvature was calculated using the height 
function (HF) technique. 

The spatial concentration profile of the converged coaxial powder flow was considered as a 
Gaussian distribution [1 1], as defined in the following equation: 

C(v , 1) — Cpg^/^([)0xp ^2 ^ 


( 2 - 1 ) 

where C is the powder concentration (the number of powder particles per unit volume), which is a 
function of the radial r and axial 1 distances in an axial- symmetrical coordinate, Cpeak is the 
powder concentration at the center of powder flow (r = 0), and Rp is the effective radius of the 
powder stream at position 1. 

The governing equations were discretized in the physical space using the finite volume method. A 
forward-staggered, fixed grid was used, in which scalar quantities were located at the geometric 
center of the cell, while velocity components were situated at the cell face centers. To discretize 
the advection terms, the flux limiter scheme MUSCL (Monotone Upstream-centered Schemes for 
Conservation Laws) [12] were applied to improve the accuracy of the upstream approximation 
and enforce the weak monotonicity in the advected quantity. The computational cycle can be 
described through the following iterative steps: 

1) Mass/momentum conservation equations and the related boundary conditions are solved 
iteratively using a two-step projection method [13] to obtain velocities and pressures. The 
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second step of the two-step projection method requires a Poisson equation to be solved for 
the pressure field: 


V 



vv 


( 2 - 2 ) 

where pn is the density from the old time step, pn+i is the pressure to be solved at the new time 

step, and ^ is the temporary velocity field computed from the first step. The density retained 

inside the divergence operator in Eq. (2-2) results in an extra term proportional to , which 

contributes to the pressure solution within the free surface transition region where ^ ^ . The 
system of linear equations formulated from the finite volume approximation of Eq. (2-2) was 
solved with an ICCG (Incomplete Cholesky Conjugate Gradient) solution technique [14]. 
Thermo-physical properties used in this step are computed from the old temperature field. 

2) The energy conservation equation is solved by a method [15] based on a finite volume 
discretization of the enthalpy formulation. Once the new temperature field is obtained, the 
thermo-physical properties are updated. 

3) Eq. (15) is solved using the PLIC-VOE [1,2,3] to obtain the updated free surface geometry of 
the melt pool and thermal field. 

4) Advance to the next time step and back to step 1 until the desired process time is reached. 

A predicted 3D track of the process is illustrated in Eigure 2.1. A melt pool is generated at the 
front of the track with the moving laser beam. The melt solidifies very quickly due to the very 
high attendant cooling rate and forms a solid track shortly after the laser moves away. The heat 
absorbed by the molten pool dissipates into the track layer and substrate quickly. The laser energy 
affects a small zone and does not degrade the material elsewhere. Corresponding cross sections 
along the scanning direction are also displayed. As can be seen, the intense heat causes substrate 
melting, which generates a good bond between the track and base material. To better examine the 
fluid motion in the melt pool and mushy zone, sequential velocity fields in cross sections along 
the scanning direction are separately shown in Eigure 2.2. The melt flows from the laser center 
area toward the periphery of the melt pool due to the Marangoni force. The fluid velocity 
deceases when the melt flows through the mushy zone where a damping effect in a porous media 
occurs. 
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(a) Cross sectional view 



(b) 3D view 



Figure 2.1 3D Deposition Profile and Temperature Distribution with Cross Section at t = 
0.11 s. Laser moves along +X direction in cross sectional and 3D view. (Beam diameter: 2.5 
mm, laser power: 416 W, travel speed: 375 mm/min, and powder mass flow rate: 6.25 

g/min.) 
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Figure 2.2 Fluid Velocity Fields in the Melt Pool in Cross Section along the Laser Scanning 
Direction at t = 0.11 s. (Beam diameter: 2.5 mm, laser power: 416 W, travel speed: 375 
mm/min, and powder mass flow rate: 6.25 g/min.) 
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3. Wire Model Introduction 


3.1. Introduction 

Based on the previous cladding code, a wire feeding model was added to this code. In the first 
step, it was necessary to construct the wire to clear the wire position. When feeding a solid wire, 
as depicted in Figure 3.1, the rectangle was assumed to be a wire. Half of the known rectangular 
width was set as r. Mq was the midpoint of the rectangular width, and was the intersection 
point of the extension of the long side and the horizontal line. The length of MqMi was set as d. 
The height to feed y and the distance to feed x also were known. As Figure 8.1 illustrates, 

01 = 02 can be achieved; therefore, according to cos0i = COS 62 , the equation ^|^ ~^y 

can be achieved. Then, d can be expressed with r, x, and y. If ^^0 < ^ ^ then point M can be 
assumed to be inside the wire position; if not, it does not belong to the wire. Finally, the wire 
position can be assumed. 



Figure 3.1 Wire position scheme 

The simulation result shown in Figure 3.2 was achieved by adding this part to the cladding 3D 
code. 
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Figure 3.2 Simulation result of wire and molten pool 


The wire shape can be seen from the simulation results. 


3.2. Analysis 

For the next step, some solver methods needed to be considered in the wire model to make the 
results more reasonable and to simulate the fluid-solid interaction. Many effects need to be 
considered when simulating this model. According to how mesh changes within a deforming 
body, three kinds of methods can be used, the Lagrangian, Eulerian, and ALE methods. 


The most important issue is to solve the solid-liquid interaction problem for the wire model,. The 
ALE (Arbitrary Lagrangian Eulerian) method can be used to handle the moving solid boundary. 

In the ALE method, the material moves in some arbitrarily specified way to enable continuous 
rezoning. Freedom in moving the computational mesh enables greater distortions of the 
continuum and more resolution than that allowed by the Lagrangian or Eulerian methods. Thus, 
the ALE formulation combines the advantages of the Lagrangian and Eulerian methods. The 
principle of an ALE code is based on the independence of the finite element mesh movement with 
respect to the material motion. 

The method aims to achieve as "Lagrange-like” behavior as possible, while simultaneously 
keeping the mesh distortions at an acceptable level. When using the ALE technique in 
engineering simulations, the computational mesh inside the domains can move arbitrarily to 
optimize the shapes of elements, while the mesh on the boundaries and interfaces of the domains 
can move along with materials to precisely track the boundaries and interfaces of a multi-material 
system. ALE-based finite element formulations can reduce to either Lagrangian-based finite 
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element formulations by equating the mesh motion to the material motion, or Eulerian-based 
finite element formulations by fixing the mesh in space. 

The sales code [16] that includes the ALE algorithm follows the general flow chart in Figure 3.3. 
The key parameter, velocity, was determined from many functions when tracing. 


celset 


♦ 

T' 

♦ 



phasel 


reerid 


Figure 3.3 General flow chart of the sales code 

Depending on the type of calculation (Lagrangian, Eulerian, or other), the shape of the mesh first 
must be modified. In the purely Lagrangian case, this allows the computing grid to follow the 
fluid motion exactly. The new set of grid vertex velocities over the entire mesh are computed 
using this equation everywhere: 


Ug Uj^ 


However, large fluid motions would create devastating effects such as contorting cells to extreme 
aspect ratios or even turning cells inside out. It is often possible to ameliorate these effects by 
moving the mesh vertices with respect to the fluid so as to maintain a reasonable mesh structure. 
For a purely Eulerian calculation, the following equation can be used: 


= 0 


Vg =0 


Whenever a vertex is moved relative to the fluid, however, the cells surrounding the vertex must 
exchange material. SALE allows a broad spectrum of rezoning possibilities by treating this 
material exchange as an advective flux. The simplest case is that of a purely Eulerian flow, in 
which the vertices are moved back to their original positions every cycle. Between this extreme 
and the Lagrangian extreme lies whatever form of continuous or discrete rezoning the user 
wishes. 
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For a continuous re-zoning, grid velocities are chosen to lie somewhere between these two 
extremes. In particular, the vertices are moved according to some relaxation rate to place vertices 
at the average position of the neighboring vertices. This usually maintains cells of reasonable size 
and proportion throughout a run. Once a set of grid velocities and Vq have been defined, it is a 
simple matter to construct the new grid and perform whatever advective flux calculations may be 
required. 

In the re grid function, the vertices are moved to new locations as specified by (U^^, V^^). Then, a 
set of relative velocities (lireh ^rei) is formed to simplify the later task of calculating advective 
fluxes. For this purpose, the vertex velocities with respect to the fluid are: 

^REL ~^G ~^L 

Then, the vertices are moved to new locations. The new grid coordinates can be computed: 

x;^^ = x"+St-u^ 

In all cases except for the purely Lagrangian calculation, the relative velocities are not zero, and 
the advective flux of mass, energy, and momentum between cells must be calculated. The flux 
calculation is performed on a cell-by-cell basis. For every cell, the volume swept out by each of 
the four faces relative to their Lagrangian positions is calculated. 

Using the concept of a cell-centered momentum flux, the average cell-centered momenta is based 
on the vertex velocities: 

umom^^ = ^ *rol^^ *(ul^^ +ul . +ul . ) 

vmom^^ = ^ *rol^^ *(vl^^ +vl . +vl . +vl,+i,,+i ) 

Then, advection changes the cell momenta as a function of the cell-centered quantities. The net 
advection changes in cell-centered momentum components are given by: 

umomp.j = f{umom,s) 
vmomp.j = f(vmom,s) 

These changes are then apportioned back to the vertices and converted to vertex velocity changes 
in the next step. Then, the vertex quantities can be updated. 

First, the new vertex masses must be calculated from the averages of the new cell masses, and a 
mass must be assigned to each vertex to obtain the time advanced velocities. In SALE, it is 
assumed that the mass in each cell is shared equally between its four comer vertices: 

M =—(M +M +M +M t 

y ^ ^ i+mj+m ^ i-mj+i/i ^ i-mj-i/i ^ i+i/ij-i/i / 
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Secondly, to adjust the velocities, the initial values set at all vertices are replaced by new time 
step values: 




m ; 


L’,j J^n+l 


Finally, the cell-centered momentum changes are distributed to each vertex of the cell: 




l^ umomp.. 


Ml 


After this section, the code returns to the beginning to start another time step using the new 
values. 

The fluid-structure interaction problems involving a free surface can be classified as follows: 

1) Rigid containers carrying a fluid with a free surface 

2) Rigid floating bodies 

3) Submerged rigid bodies moving close to a free surface 

When modeling the wire feed process, the first step is to find the equations of motion for the 
structure as functions of the interaction pressure. The second step is building a boundary- value 
problem for the pressure as a function of the acceleration of the structure. The third step deals 
with coupling the pressure solution with the equations of motion of the structure. 


When implementing the ALE method in the wire feed, the general method was simplified. For the 
case of fluid-rigid body interaction, the wire was modeled as a rigid body. For it to be classified 
as merged solid and non-merged solid, in the proposed wire model, the nonmerged solid and free 
surface modeling need to be taken into account. A numerical method that can simulate this 
interaction requires coupling between the fluid solver and a code for the structural dynamics. 

Just like the original code, the fluid is modeled with the Navier-Stokes equations. The mass and 
momentum conservation equations are as follows: 

VV = 0 

^(pV) + V -(piV -%)V) = V -V^) + pg + S, 

dt p, K p, 

The free-surface dynamics are tracked using the Volume of Fluid (VOF) method. The VOF-type 

interface-capturing method is used to compute the deformation of the free surface: 

dF ^ ^ 

— + (y-L)-VF=0 
dt 

The motion of the rigid body in the six degrees of freedom is described by the equations of linear 
and angular momentum, set in the inertial reference frame, and given by the following two 
equations: 
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MX =F 

RIR-'a + hxRIR-^a = Mg 

Here, M is the boat mass, ^ is the linear acceleration of the center of mass, P represents the 

force acting on the boat, ^ and ^ are the angular acceleration and velocity, respectively, 
is the moment with respect to ^ acting on the boat, I is the tensor of inertia of the boat about 

the body-fixed reference system axes, and R is the transformation matrix between the body-fixed 
and the inertial reference systems. 

The equation of variation of angular momentum in the form referred to as the center of gravity is: 

where ^ and ^ are the absolute angular acceleration and angular velocity, respectively, is 

the total moment with respect to ^ , and is the tensor of inertia of the body about the 

(G,x, y, z) axes. 


To simulate the motion of bodies floating freely at the free surface, the equations of motion of the 
rigid body can be solved together with the flow solver. In the coupling with the flow solver, the 
6DoF dynamical system receives, at each time step, the value of the forces and moments acting 
on the rigid body, and it returns the values of the new position, as well as linear and angular 
velocity. In the flow solver, these data are used to update the computational grid by a mesh 
motion strategy based on elastic analogy. The flow equations in ALE form are solved in the new 
domain. 

The general idea is to decouple the problem in the following way. The flow solver computes the 
flow around the body as usual. Taking into account the fluid viscosity, flow turbulence, and 
deformation of the free surface, the forces and moments acting on the body are then calculated by 
integrating the normal (pressure) and tangential (friction) stresses over the body surface. 
Following this, the rigid body equations can be solved using these forces and moments as input 
data at each time step, and the motion accelerations, velocities, and displacements (translations 
and rotations) are obtained by integrating in time. The position of the body is then updated, and 
the fluid flow is computed again for the new position. In the flow solver, these data are used to 
update the computational grid, and the flow equations in ALE form are solved in the new domain. 
By iterating this procedure over time, the body trajectory is obtained. 

However, many problems still need to be solved. For example, the ALE differential form of the 
conservation equations for mass, momentum, and energy need to be rewritten in another form. 
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All one must do to obtain the ALE form of the conservation equations is to substitute in the 
various convective terms, such as the material velocity with the convective velocity. The ALE 
differential and integral forms of the conservation equations will be reset and used as a basis for 
the spatial discretization of problems in the fluid dynamics and solid mechanics. 

In addition, the computer implementation of the ALE technique requires the formulation of a 
mesh-update procedure that assigns mesh-node velocities or displacements at each time step of a 
calculation. The remesh algorithm strongly influences the success of the ALE technique. 

Also, mesh regularization is necessary to keep the computational mesh as regular as possible 
during the entire calculation in order to avoid excessive distortions, squeezing of the computing 
zones, and to prevent mesh entanglement to decrease the numerical errors caused by mesh 
distortion. 

The boundary conditions are of utmost concern in the discussion of the specificities of the ALE 
formulation in fluid dynamics. The same boundary conditions employed in Eulerian or 
Lagrangian descriptions are implemented in the ALE formulation. Along the boundary of the 
domain, kinematical and dynamical conditions must be defined. Along solid-wall boundaries, the 
particle velocity is coupled to the rigid structure. Due to the coupling between the fluid and 
structure, specific conditions are needed to ensure that the fluid and structural domains will not 
detach or overlap during the motion. These coupling conditions depend on the fluid. Eor a viscous 
fluid, the coupling between the fluid and structure requires that velocities coincide along the 
interface. In practice, two nodes, one fluid node and one structural node, are placed at each point 
of the interface. The fluid is treated in the ALE formulation; the movement of the fluid mesh may 
be chosen completely independently of the movement of the fluid itself. The equations for a 
Newtonian fluid in terms of stress and the strain rate need to be solved, and in three dimensions, 
the strain rate in each cell should be calculated from the velocity gradient, taking other equations 
into account. 

3.3. Future Work 

The condition where the wire is implemented just above the melt pool, instead of feeding a solid 
wire into contact with the liquid based on previous cladding code, will need to be simulated. The 
complete wire feeding model needs to be done in the future work. 
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4. Parallel Computing 


A Message Passing Interface (MPI) was designed for high performance on both massively 
parallel machines and workstation clusters. MPI, now being used for parallel computing, is a 
specification for the developers and users of message-passing libraries. The goal of MPI is to 
provide a widely used standard for writing message-passing programs. The interface attempts to 
be practical, portable, efficient, and flexible. MPI specifications have been defined for C/C++ and 
Fortran programs [17]. 

To parallelize the ICCG solver, using the parallelize version of the PCR solver [18] as a guide, 
the ICCG code was divided and then grouped into the following 14 different files: surflO.cpp, 
mpi_slice_moves.cpp, init_funcs.cpp, iccg.cpp, functions.h, functions.cpp, function_prototypes.h, 
extern_globals.h, convection.cpp, constants.h, comp_cal.cpp, claddingSDmpi.cpp, bc_funcs.cpp, 
makefile. Each file consists of related function(s)/variables. For instance, the integer variables are 
saved in the “constants.h” file, while float variables are saved in the “extem_global.h” file. The 
3D arrays F, FN, and FP were identified as those that require two slices, so they were indexed 
with the “fbi” variable. Similarly, the ID arrays DX, RXD, and XI were indexed with the “gi” 
variable, indicating that each node must have its own copy of each of these arrays. 

In cases in which the for() loop was once for (i=0; i<IMAX; i++), the loop has been changed to 
for (i=imaxBeg; i<imaxEnd;i-i-i-)'^^^l Similarly, in cases in which the for() loop was once for (i=l; 
i<IMl; i++), the loop has been changed to for (i=imlBeg; i<imlEnd;i-i-i-)^^^l 
The complexity of indexing in the ICCG solver also was reduced by categorizing nodes into three 
different natural types: node_0, mid-nodes, and the end node. 

if(myRank = = 0) f 

nodePos = NODE_0; 
imaxBeg=0; imaxEnd=sliceN; 
imlBeg-1; imlEnd=sliceN; 

} 

else if (myRank< numMPINodes-1) { 
nodePos = NODE_i; 
imaxBeg=l; imaxEnd=sliceN+ 1 ; 
imlBeg=l; iml End=sliceN+ 1 ; 

} 

else if (my Rank = = numMPINodes-1) { 
nodePos = NODE_N; 
imaxBeg-1; imaxEnd=sliceN+ 1 ; 
imlBeg-1; imlEnd=sliceN; 
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Finally, at the end of each for() loop, the sliceMove() function was used to update the ghost slices 
for those arrays whose values had changed. 

The results from the temperature files retrieved from the current paralleled version are as follows: 



Figure 4.1 Results from the serial code 



Figure 4.2 Picture of dividing model for the serial code 
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Figure 4.3 Results of the left half model for the parallel code 



Figure 4.4 Results of the right half model for the parallel code 


The meshing grid pictures show that there was one layer above the hottest grid that was within 
one grid. The interpolation method was used to calculate the temperature across the melt pool, so 
the extra layer was derived using the mean value of room temperature (blue part in the figures) 
and the highest temperature (red part in the figures). 

The code appears as follows: 

F_SUM3 = F_SUM3 + F[fbi+3][M][N]; 

F_SUM3 = F_SUM3 + F[fbi-3][M][N]; 
if(FP[fbi+3][M][N]>=EMFl && F[fbi+3][j][k] < EMF) 

FOFPPP = 0.0; 
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if(FP[fbi-3][M][N]>=EMFl && F[fbi-3][j][k] < EMF) 

FOFMMM = 0.0; 

The variables “F” and “FP” both require three slices. In the allocation section, TotSlicesS was 
added for processes that require one or two ghost slices. 

if(nodePos == NODE_0 II nodePos == NODE_N) 

{ 

// end procs need 1 ghost slice 
TotSlices-sliceN+1; 

TotSlicesl =sliceN+2; 

TotSlices3 = sliceN+3; 

} 

else 

{ 

// all others need to keep 2 ghost slices 
TotSlices-sliceN+2; 

TotSlices2 =sliceN+4; 

TotSlices3 =sliceN+6; 

} 

For cases in which nodePos does not equal NODE_0, fbi=i+l has been changed to fbi=i+2, as 
follows: 

if (nodePos == NODE_0 ) 

{ 

fbi=i; gi=i+sliceO_G; 

} 

else 

{ 

fbi=i+2; gi=i+sliceO_G-l; 

} 

The MPI functions were compiled on the cluster successfully. Some C++ needs to be included in 
directives to the header of the makefile to complete the successful compilation of the code. 
However, some “subscript out of range” issues still remain that cause the solver to crash. These 
issues existed in the sequential version of the ICCG solver and transferred to the parallel version. 
Some functions need to be synchronized across the different cores. We wrote the following 
function to achieve this functionality: 

void synchronizeFlag(int &flag, int syncStyle) { 

MPI_Barrier(MPI_COMM_WORLD ); 
bool IsSame = true; 

int HempFlag - (int ^)malloc(sizeof(int)^l); 

tempFlag[0] =flag; 

int ^flagsVector = NULL; 

if(myRank -- 0) { 

flagsVector = (int ^)malloc(sizeof(int) ^ numMPINodes);} 


20 



MPI_Barrier(MPI_COMM_WORLD ); 

MPI_Gather(&tempFlag[0], 1, MPIJNT, flagsVector, 1, MPIJNT, 0, 
MPI_COMM_WORLD ); 

MPI_Barrier(MPI_COMM_WORLD ); 

ifimyRank = = 0)f 

for(int i=l; i<numMPINodes; i++){ 
if(flagsVector[0] !=flagsVector[i]) 

IsSame = false;} 
ifilsSame !- true){ 
switch( syncStyle ){ 
case 0: 

tempFlag[0] = computeMinimumFlag(flagsVector); 
break; 
case 1: 

tempFlag[0] = computeMaximumFlag(flagsVector); 
break;}}} 

MPI_Barrier(MPI_COMM_WORLD ); 

MPI_Bcast(tempFlag, 1, MPIJNT, 0, MPI_COMM_WORLD); 
flag = tempFlag[0]; 

//free stuff 
free( tempFlag ); 
free(flagsVector); 

} 

Using the above function, the following flags were synchronized: FLGMELT, FLGVLCT, 
DELT, FLGC, NFLGC, and NOCON. 

We tested heat conduction in both the serial and parallel codes under the following parameters: 

Table 4.1 Simulating parameters 


Parameter 

Value 

Domain Length (mm) 

3.2 

Domain Height (mm) 

1.84 

Domain Width (mm) 

3.2 

Substrate Height (mm) 

1 .47 (80% of Domain Height) 

Element Size (pm) 

40 

Laser Power (watts) 

1000 

Pulse Width (secs) 

0.001 
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The results were as follows: 



(a) Result from serial code 



Figure 4.5 Results after 0.0013secs 



(a) Result from serial code 



(b) Result from one core of the parallel code 


Figure 4.6 Results after 0.021secs 


22 



(a) Result from serial code (b) Result from one core of the parallel code 

Figure 4.7 Results after 0.0013secs 



Figure 4.8 Cooling curve (serial and parallel results overlapped) 


Figure 4.5 through Figure 4.8 indicate that the results from both the serial and parallel codes were 
similar. 
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5. Smoothed Particle Hydrodynamics Model 


5.1. Abstract 

Smoothed Particle Hydrodynamics (SPH) is a new numerical tool useful for hydrodynamic 
simulations. Compared to traditional techniques such as the VOF method, SPH is a more 
convenient and robust technique for complicated simulations involving fluids. In this report, the 
modeled SPH is described briefly, and its advantages are noted through a comparison with the 
VOF model. The working principles of SPH are discussed, and the implementation of this 
technique in the python code is discussed in detail. The step-by-step flow of the program is 
discussed, and results are plotted. 

5.2. Theory 

The primary advantage of SPH over other computation techniques is that it is a mesh-free 
Lagrangian method. In the Lagrangian description, the position and physical properties of 
particles are described in terms of the material or referential coordinates and time. The SPH 
method works by dividing the fluid into a set of discrete elements, referred to as a particle. These 
particles have a spatial distance (known as the "smoothing length" and typically represented in 
equations by h) over which their properties are "smoothed" by a kernel function. This means that 
the physical quantity of any particle can be obtained by summing the relevant properties of all the 
particles that lie within the range of the kernel [19]. 
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Figure 5.1 Pictorial representation of the kernel function and smoothing length 

The contributions of each particle to a property are weighted according to their distance from the 
particle of interest and their density [20]. Mathematically, this is governed by the kernel function 
(symbol W). Commonly used kernel functions include the Gaussian function and the cubic spline. 
The latter function is exactly zero for particles further away than two smoothing lengths (unlike 
the Gaussian, for which there is a small contribution at any finite distance). The latter function 
has the advantage of saving computational effort by not including the relatively minor 
contributions from distant particles. 
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Figure 5.3 Cubic spline curve 


To utilize the full power of SPH, the smoothing length can be varied over time. For example, the 
current simulation has two smoothing length values, one for relatively short length ranges, which 
computes forces such as compression, friction, etc., and one for long ranges, which computes 
forces caused by surface tension, etc. In the current simulation, only the fluid particles are 
considered as SPH particles for computational efficiency [21]. 



Figure 5.4 3D representation of the kernel function 
5.3. Comparison with VOF Technique 



Figure 5.5 Non-fluid and fluid interface in VOF technique 


The Volume of Fluid (VOF) technique uses fluid fractions to determine the fluid cells and solid 
cells. This technique primarily uses the Eulerian approach (fixed-grid method). The main 
disadvantage of the Eulerian approach is that it is difficult to conserve mass, as the cells have no 
mass of their own. However, in the SPH technique, each cell possesses its own mass, which is 
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conserved throughout the system. In the VOF technique, the motion of the fluid is computed from 
the momentum equations, including the Navier-Stokes equation, and many conditions to maintain 
continuity. The SPH technique eliminates all of these complicated equations, replacing them with 
much simpler force and velocity calculations [22]. Calculating the advection and solving the 
Pressure-Poisson equation take most of the computation time and make real-time processing 
unachievable. Separate solvers must be implemented to solve the inverse matrix calculation for 
the pressure matrix. In the VOF technique, the transport phenomenon by advection is implicit, 
requiring nested loops in three dimensions and taking most of the computation time. The VOF 
technique is very difficult to parallelize, as processors of the machine will be waiting for values 
from the other processors to finish computing. The SPH model is highly parallelizable because 
the fluid is considered as particles with their own physical properties. The primary advantage of 
the SPH model is that most of the calculation can be done explicitly. Hence, array operations, 
which are faster, can be used instead of nested loops. 
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5.4. Program Flow 



Initialize the global variables and arrays used for 
computation 



Create substrate and solid objects in the domain 


Create arrays for the substrate with intiial condition, 
and create wire 


W 


Material properties and initial conditions 


Import the material properties of the required 
material and pass it to the variables 



Write results to the file 
after a predefined 
number of iterations 


Update enthalpy 
based on the heat 
obtained 



Apply heat 
condcution, radiation 
losses and update 
enthalpy 


Process all 
movement in the 
simulation 


Calculate 

temperature from 
the enthalpy 


Calculate solid-fluid 
interaction forces 


/ 


Check melting and do 
phase change ' 
calculations 


Calculate fluid forces 
based on SPH 
technique 
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6 . SPH Model Implementation 


6.1. Initialization 

The global variables are set to their initial conditions in the beginning of the code. The dimension 
of the domain is specified with a desired resolution from which the grid size is calculated in three 
dimensions. The time step at which the iteration is to be calculated is specified and declared as a 
global variable, which adds to the total simulation time in each cycle. The arrays used are 
“ctypes”, which have the efficiency of C arrays and can be used in Python to do calculations 
while reducing memory usage. 

The initial position of the laser in the domain, the velocity and power of the laser, and the type of 
distribution of power are specified as part of the initialization. Two types of laser power 
distributions exist in the current code: “top hat” and “Gaussian” beam. Figure 6.1 shows the 
distribution of each laser beam. 



Top-hat distribution Gaussian distribution 

Figure 6.1 Two types of laser beam distribution 


A laser map that has the indices of each column of the laser beam in the domain is created, with 
the power in that column according to the profile. This yields a more realistic “ray traced laser.” 
With the ray traced laser, the solid objects hit by the laser as it moves downward to the substrate 
are taken into account. This projected laser will cast a shadow of the solid that obstructs the rays 
in the substrate, which is more realistic than applying laser power onto the substrate alone. This 
approach is better explained in Figure 6.2, which shows the laser projected onto a substrate with a 
solid in between, with and without ray tracing. 
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Without ray-tracing With ray-tracing 

Figure 6.2 Laser projection on substrate scheme 

6.2. Creating Solids and Applying Material Properties 

The material properties are stored in a dictionary format for each material. This will allow more 
materials to be added into the file without changing the original code. The material used in the 
current code is Ti6A14V and its properties. The next step is to create the solids used in the 
simulation. First, a substrate having the desired dimensions is created and placed in the domain 
with a user-defined offset. The program has functions that create cylinders, cuboids, spheres, 
frustums, discs, and other shapes. This shows that a solid of any shape can be added to the 
simulation. The substrate is declared as a fixed, solid object so that it will not be moved according 
to forces such as gravity. The material properties are applied to the solid, and the substrate is 
placed under ambient temperature. A temperature-enthalpy lookup table is created as a part of the 
initialization and used in the temperature calculation. The mass of the object is calculated based 
on the mass of each cell, which is calculated from the solid density and volume. 

Wire is created in the same way substrate is created, as explained previously. The additional input 
values needed to create the wire include the diameter and length of the wire, which should not 
exceed the size of the domain. The next input value for the wire is its angle with respect to the X- 
Z plane and the X-Y plane. The velocity of the wire is a user-defined value, and the gravity of the 
wire is restricted to prevent the wire from falling on the substrate due to its own weight. 

6.3. Iterative Calculations 

The iterative calculation done at each time step is divided into 10 different steps. 
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6.3.1. Laser Projection 

As explained previously, a ray traced laser is projected to the domain in the Z-axis. The heat from 
the laser is added to the enthalpy array at the intersection positions. The intersection position is 
obtained by tracing each ray down its column to find the first solid cell in its path. 

6.3.2. Heat Conduction 

The conduction calculation begins with the creation of a conductivity matrix. The conductivity 
value is given as a material property based on the temperature. The conductivity matrix is 
calculated through the interpolation of the range of values of “K” with the temperature. The next 
step is to create a gradient matrix based on the differences in temperature in three different axes. 
The heat due to conduction is then updated on the enthalpy matrix. 

6.3.3. Radiation Losses 

Radiation losses are subtracted from the skin cells of the solids because radiation is a surface 
phenomenon. To find skin cells, scipy-ndimage-morphology libraries are used, which include 
functions such as binary erosion and binary dilation for array operations. The enthalpy array is 
updated after radiation losses. Both the conduction and radiation loss function are performed in 
parallel, as the calculation of one does not depend on the other. This increases the speed and 
efficiency. 

6.3.4. Temperature Calculation 

The temperature is calculated from the enthalpy array, which was modified by the laser power, 
conduction, and radiation. The temperature array is computed by interpolating the enthalpy values 
with the temperature in the previously calculated temperature-enthalpy lookup table. Using the 
lookup table makes the computation much faster and avoids values exceeding a reasonable limit 
with larger time steps. 

6.3.5. Phase-Change Calculation 

The phase-change calculation is divided into three parts: checking melting, solidification, and 
creating voxels. 

a) Checking melting 

The melting process begins with checking the temperature array and finding the indices of cells 
exceeding the melting temperature of that material. The indices are mapped with a temporary 
array to form new fluid cells. The function that creates new fluid cells creates a list of fluid 
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positions, fluid indices, fluid velocities, and fluid forces. The fluid velocities are inherited from 
the solid cells, and forces are declared as zero when created because they will obtain their value 
in the next step. 

b) Solidification 

The liquid cells can solidify in the simulation due to conduction, radiation, and movement of the 
laser. The fluid cells are mapped with the temperature array to check for cells with temperatures 
less than the solidus temperature of the material. These cells are deleted from the fluid lists and 
added back to the solid cells. 

c) Creating voxels 

The final step in a phase change is to create the arrays required for fluid computation based on the 
results given by melting and solidification. The values computed by this function are used later to 
create SPH particles and fluid movement. Arrays such as fluid space are created and the fluid 
positions determine the location of the cells. This will indirectly help the movement of the fluid 
because the fluid position is the variable that will be updated in the move function, which is 
explained later. 

6.3.6. Updating Smoothing Length 

The smoothing length is calculated here from a kernel function, which is cubic spline in the 
current model because of computational efficiency. This step in the calculation can be 
considerably slower as more fluid cells are added into the simulation. Hence, parallel computing 
requirements are high in this step. For efficient parallelization, a list consisting of reference 
numbers of fluid particles and their positions is created. This list will have multiple lists of fluid 
particles, considering two particles at a time. A cubic spline worker function calculates the 
smoothing length for two interacting particles. The list is passed to a multiprocessing pool that 
conducts the above operation in parallel, saving computation time and using resources efficiently. 
The final output of this function is a list of lists, with each list having the position of the two 
interacting particles, the direction of interaction, smoothing length 1 for a short range of forces, 
and smoothing length 2 for a higher range of forces. 

6.3.7. Calculating Fluid Forces 

In this step, the SPH model is applied to the fluid particles. The fluid forces, which were declared 
as zero in the phase change step, are changed to gravity forces in this step. The list created in the 
previous step is used in this function. Each pair of particles is taken one at a time to find forces 
such as the compression force, friction force, and force due to surface tension. Smoothing length 
1 is used for compression and friction forces, and smoothing length 2 is used for the surface 
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tension force. Material properties such as stiffness, viscosity, and surface tension are taken into 
account in finding these forces. The forces are added to one fluid particle and subtracted from the 
other based on the direction of the force. 

6.3.8. Solid-Fluid Interaction Calculation 

In this step, the interaction forces between the solid and the fluid are computed. The first step in 
this calculation is to find the neighboring cells of each fluid particle in three dimensions and 
check if they are solid cells. The smoothing lengths are calculated for the interacting pairs in the 
same way they were calculated in the previous steps. The reaction forces and friction forces are 
computed using the two smoothing lengths, respectively, and added in the direction of action. The 
forces added in the fluids are subtracted from the solid body. 

6.3.9. Movement 

The movement function can be divided into four different motion calculations: 

a) Moving solids 

In this step, each solid body is moved with its predefined velocity. In the current simulation, the 
substrate is stationary, and the wire moves with a velocity in the angle in which it was created in 
the domain. Movement due to gravity forces is computed in this step, but because the gravity is 
restricted for the wire and the substrate is fixed, the force of gravity has no influence on this solid. 
However, the gravity force becomes a critical factor in fluid and simulating powder problems. 

b) Collision 

The next step in the movement calculation is to check for collision. Each solid has a parameter 
named “target” that determines whether the solid is to be considered for collision calculation. 

This condition is checked in the beginning and, if true, is passed to a function that checks for 
collision by determining if the space arrays of two solids overlap; any overlapping indices are 
returned for further calculations. In the current simulation, if two solids collide (i.e.; the wire and 
the substrate) then the wire projection is incorrect or the wire is moving so fast that it hits the 
bottom of the melt pool, which in the real world is considered a failure. The code breaks at this 
point, giving an error message that there has been a collision. However, further calculations are 
written into the code for powder simulation, which ideally involves collision. In such cases, if a 
particle does collide with a solid, the moving solid is backed up by half a time step, and the 
velocity is recalculated based on the coefficient of restitution of the material and averaging the 
mass. This calculation is used for a more realistic solid-to-solid collision, which causes 
movements such as bouncing and others. 

c) Fluid movement 
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The third step is to process the motion of the fluid particles in the simulation. The fluid velocities 
are updated with the forces calculated with SPH smoothing. The fluid positions and indices are 
updated from the velocities. The movement part will be updated automatically in the create voxel 
junction based on the updated positions, 
d) Laser movement 

Finally, the position of the laser beam is updated based on the user-defined speed of the laser and 
the time step. 

6.3.10. Increment Time and Save Data 

The final step in the iteration is to increment the simulation time by one time step and to repeat all 
of the above steps until the simulation time reaches the user-defined time. The results, including 
the simulation time, temperature, and phase arrays are saved into a data file at a predefined 
number of iterations. These data files are used to create the images that are assembled to form the 
simulation movie. 

6.4. Results 

Results from the simulation are plotted using visualization software such as Mayavi and Matplot. 
Tests are conducted for different input parameters, and the results are plotted. 

1) Resolution 

Results are plotted with varying resolution values by keeping all other parameters constant. The 
results appear in Figure 6.3. 


0.000020 s I 0,000020 s 


a) Resolution - 0.0001 b) Resolution - 0.00006 
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0,000020 s 


c) Resolution - 0.00003 

Figure 6.3 Result under various resolutions 

As shown in Figure 8.3, the lower the resolution, the more cells and the finer the mesh. However, 
lower resolution increases the computation time. 

2) Wire Diameter 

Results are plotted with varying wire diameters, as shown in Figure 8.4. The domain size is 
1cm* 1cm* 1cm, with a resolution of 0.00006 


O.Q00020 s 

0. 000020 s 

a) Diameter - 1mm 

b) Diameter - 3mm 

Figure 6.4 Result using various wire diameters 


3) Wire Length 

Images of varying wire lengths are shown in Figure 8.5. The domain size and resolution are the 
same as in the previous result. 
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0.000020 s 


0,000020 s 


a) Length - 3mm b) Length - 5mm 

Figure 6.5 Result using various wire lengths 


4) Wire Angles 

The wire angle can be varied in the X-Z plane and the X-Y plane. Results are plotted with 
different values of X-Z and X-Y angles in degrees. As shown in Figure 8.6, the wire can be 
rotated about a plane in a specified angle. There is separate input options for the X-Z plane and 
the X-Y plane angle of rotation. 


0.000020 s 


a) X-Z - 70, X-Y - 0 


0,000020 s 


b)X-Z- 70, X-Y - 30 
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0.000020 s 


c) X-Z- 70, X-Y - 60 

Figure 6.6 Result using various wire angles 

5) Wire Offset 

The offset is the dimension in the domain at which the wire is placed after it is created with the 
desired length, diameter, and angles. Results are plotted with different offset values. The offset 
from the domain at which the wire has to be placed can be specified. As shown in Figure 8.7 the 
offsets specified makes the wire appear in different location when the simulation starts. The 
figure is shows a 10mm, 10mm, 10mm size domain. 


0.000020 s 

0.000020 s 

a) Offset - 3mm, 5mm, 6mm 

b) Offset - 5mm, 6mm, 8mm 

Figure 6.7 Result using various wire offsets 

6.4. 1 . Temperature Plots Using Matplot 



The temperature array is plotted using Matplot, as shown in 8.8. 
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Figure 6.8 Temperature(K) plot using Matplot 

In the figure, the X and Y axes represent the number of cells in the domain, and the color bar 
shows the temperature range and the color. The laser started from the center of the substrate and 
reached halfway to the end of the wire. The figure clearly shows the shadow of the wire falling on 
the substrate due to ray tracing. The wire is hotter than the substrate due to having less cross- 
sectional area. 


6.4.2. Phase Plots Using Mayavi 


0.009300 s 


Figure 6.9 Phase plots showing melting 
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Figure 6.9 shows the phase plots with the melt pool formed. The title is the time of the simulation 
at which the image is plotted, 0.009300(sec). As shown in the figure, the tip of the wire is starting 
to melt and drip off from the wire to the melt pool. 

6.5. Conclusion 

Smoothed particle hydrodynamics has been shown to be one of the most efficient ways to 
numerically simulate problems involving multiphase flow and phase changes. Its computational 
efficiency and simple architecture makes it more robust and flexible than other simulations. The 
current simulation is a partial SPH model that considers fluid as SPH particles. A full SPH model 
would allow for more realistic calculations, such as the bending of the wire and other parameters. 
The technique could be integrated with CAD software in the future to easily create the solids and 
depositing parts of complicated geometries. 
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7. Numerical Analysis of Thermal Stress and Deformation in 
Multilayer Laser Metal Deposition Processes 


7.1. Introduction 

Highly localized heating and cooling during the Direct Metal Deposition (DMD) process 
produces nonuniform thermal expansion and contraction, resulting in a complicated distribution 
of residual stresses in the heat affected zone and unexpected distortion across the entire structure. 
The residual stresses may promote fractures and fatigue and induce unpredictable buckling during 
the service of deposited parts. The distortion is often detrimental to the dimensional accuracy of 
the structure; therefore, it is vital to predict the behavior of materials after the DMD process and 
optimize the design/manufacturing parameters to control the residual stresses and distortion. 

The focus of this research was to investigate both the thermal and mechanical behavior of the 
DMD process of Stainless Steel 304. Based on the finite element analysis package ABAQUS, a 
3D, sequentially coupled, thermo-mechanical model was developed to simulate the transient 
temperature field, residual stress, and final deformation. A laser displacement sensor was used to 
record the deflection of the substrate caused by thermal stresses during the deposition process. By 
comparing the experimental results with the simulation results, the numerical model was 
validated. 

7.2. Thermal Analysis 

The numerical modeling consisted of two main steps. A transient thermal analysis first was 
carried out to generate the temperature history of the entire work piece, which then was used as 
input in the subsequent mechanical analysis. 

7.2.1. Initial and Boundary Conditions, and Energy Distribution 

The initial condition applied to the model was the ambient temperature, which was set as room 
temperature, 298.15 K. The boundary conditions, including thermal convection and radiation, 
were considered. A circular laser beam is shot onto the substrate vertically with a constant and 
uniform power density. Thus, the heat source team Q was considered a constant and uniformly 
distributed surface heat flux, defined as follows: 
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where a is the absorption coefficient, P is the power of the continuous laser, and r is the radius of 
the laser beam, a was set as 0.4 according to numerous experimental results, and r was set as 1.25 
mm. 


7.2.2. Powder Addition 

During modeling, the continuous powder addition process was divided into many small time 
steps. In each time step, using the “Model Change,” a set of elements was added onto the 
substrate to form rectangular deposits along the centerline of the substrate (Figure 7.1). The width 
of the deposits was assumed to be the same as the diameter of the laser beam, and the thickness of 
the cladding was calculated from the laser travel speed and powder feed rate with an efficiency of 
0.3. The geometry of the deposits was updated at the end of each step to simulate the 
corresponding boundary conditions. 


7.2.3. FEA Model 


As shown in Figure 7.1, a finite element model for a 1-pass, 3-layer DMD process was built. Two 
cases were simulated with different process parameters, including laser power, laser travel speed, 
and powder feed rate, as shown in Table 7.1. These parameters were chosen according to the 
criterion that the final geometry of the deposits and the total energy absorbed by the specimen be 
the same in both cases. 




Figure 7.1 Dimensions of DMD specimen. The unit is millimeter. 
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Table 7.1 DMD process parameters 


Case number 

Laser power (fk) 

Laser travel speed {mm/ min) 

Powder feed rate (g/min) 

1 

607 

250 

6.3 

2 

910 

375 

9.4 



Figure 7.2 Meshing scheme 


7.2.4. Numerical Temperature Field 

Figure 7.7.3 shows the temperature field of the melt pool and the surrounding areas at different 
times in Case 1. The peak temperature during the process was 2500 K, while the simultaneous 
lowest temperature was close to room temperature. Considering the large temperature differences 
and small geometrical dimensions, very large temperature gradients would exist in the DMD 
process. 
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Figure 7.7.4 Contour plots of temperature field at different times (Case 1) 


Figure 7.5 shows the temperature of the nodes in the x (length) and y (height) directions in Case 1 
at t = 2.7 s. The x-direction nodes were selected on the top surface of the substrate, while the y- 
direction nodes were selected along the height of the deposits. The top surface of the substrate 
had a maximum temperature of 1 1 10 K right below the center of the laser beam, and the 
temperature decreased gradually along the x direction. In the y direction, the deposits had a 
maximum temperature of 2248 K on the top surface, and the temperature decreased rapidly to 
1 1 10 K. The slopes of the temperature curves clearly represent the thermal gradients in the x and 
y directions. In the x direction, the temperature gradient had a maximum value of 450 K/mm; in 
the y direction, the maximum temperature gradient occurred near the top surface of the deposits 
and was as large as 1400 K/mm before decreasing along the negative y direction. These steep 
thermal gradients would induce large compressive strains within the deposits and substrates. 
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7.3. Residual Stress 


The residual stress distribution within the final deposits appears in 
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(c) (d) von Miscs Stress 

Figure 7.6 (one-quarter of the deposits are hidden to show the internal residual stress). Normal 

stresses On, O22, <^33 along three spatial directions appear in 
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Figure 7.6(c), respectively, and the von Mises stress is shown in 
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(C) (d) von Mises Stress 

Figure 7.6(d). The figures show that compressive residual stress existed in the top free surface of 

the deposits, which was caused by the steep temperature gradient. As Figure 7.5 indicates, the 

temperature gradient in the x and y directions reached as high as 450 K/mm and 1400 K/mm, 

respectively. The expansion of the hotter top layer was inhibited by underlying material, thus 

introducing compressive stress in the top surface. Residual stresses in the lower part of the 

deposits were mostly tensile stresses due to the cool-down phase of the molten layers. After the 

deposition was finished and the laser was turned off, the remelted lower part of the deposits 

began to shrink, and this shrinkage was restricted by the underling material, thus inducing tensile 

stresses. 
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Figure 7.6 Contour plots of residual stress field within deposits 


The distribution and magnitude of the residual stresses in the top surface of the deposits appear in 
Figure 7.7. In the x direction, the middle part of the top surface was compressed with a stress 
magnitude of approximately 200 MPa, while the two edges along the z direction were slightly 
tensioned. In the y direction, the residual stresses almost vanished. For the normal stresses along 
the z direction, tensile stresses with magnitudes of approximately 100 MPa existed near the 
center, and compressive stresses ranging from 0 to 200 MPa existed near both ends. 



(a) al 1 in top surface (b) o22 in top surface 
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Figure 7.7 Residual stress in top surface of deposits 


7.4. Deformation 

During the DMD process, the substrate will continuously experience expansion and shrinkage, 
ultimately maintaining a deformed shape. In this research, deflection in the y direction was the 
main deformation under consideration and was observed in both the experiments and simulations. 

7.4.1. Experimental Setup 

As shown in Figure 7.8, in the experiment, the substrate was clamped at the left end to prevent 
rigid body motion. A Keyence’s LK-G5000 series laser displacement sensor with the accuracy of 
0.25 micro shown in Figure 7.9, was placed directly below the right end of the substrate to record 
the displacement of the free end in the y direction with a frequency of 25 Hz during the process. 
The experimental results appear in Figure 7.10. 




Dtsplocem^nt 

WMSurement 


Figure 7.8 Experimental setup for measuring displacement during deposition 
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Figure 7.9 Laser displacement sensor 

Figure 7.10 shows the comparison of the deflection of the substrate between the experimental and 
simulation results for both cases. These plots indicate that the trend of the deflection calculated 
from the simulation matched very well with the experimental results. For each deposition layer, 
the substrate first bent down due to thermal expansion on the top surface and then bent up due to 
thermal shrinkage during the cooling process. After completely cooling down, the substrate 
maintained the deformed shape. 
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(a) Deflection in Case 1 
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(b) Deflection in Case 2 

Figure 7.10 Simulation and experimental results of deflection of substrate 


The differences in the final deflection value between the simulation and the experiment were 
28.5% and 24.6% for Case 1 and 2, respectively. Several possible reasons could exist for these 
differences. First, errors existed in experiment setup. In the simulation, the laser beam traveled 
exactly along the centerline of the substrate. However, this could not be accomplished perfectly in 
the experiments. These offsets would affect the deflection to a large extent because the deflection 
is sensitive to the position of the heated zone and measuring point. Secondly, the laser 
displacement sensor does not track the displacement of one particular node. Rather, it works by 
sensing the signal reflected by the obstacle, so the positions it monitors are always changing as 
the substrate continues to deform. The simplifications and assumptions considered in both the 
thermal and mechanical analysis are also important factors contributing to the differences. 


7.5. Conclusion 

A sequentially coupled, thermo-mechanical, finite element model was developed and validated 
experimentally. The results showed the temperature distribution, residual stress, and deformation 
within the formed deposits and substrates. Finite element modeling can be used to effectively 
predict the resulting mechanical behavior of materials after Laser Aided Direct Metal Deposition 
processes. More efforts are needed to explore more complicated situations in industry and to 
optimize design/manufacturing parameters to control the residual stress and distortion. 
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8. Solidification Microstructure 


8.1. Cellular Automaton (CA) - Finite Element (FE) Model Introduction 

A predictive model based on a Cellular Automaton (CA) - Finite Element (FE) method [23] was 
developed to simulate microstructure evolution during metal solidification for a laser-based 
additive manufacturing process. The macroscopic FE calculation was designed to update the 
temperature field and simulate a high cooling rate. In the microscopic CA model, heterogeneous 
nucleation sites, preferential growth orientation, and dendritic grain growth kinetics were 
simulated. The CA model was able to show the entrapment of neighboring cells and the 
relationship between undercooling and the grain growth rate. The model predicted the dendritic 
grain size and morphological evolution during the solidification phase of the deposition process. 
The model parameters for the simulations were based on stainless steel 316(SS316) and T16A14V. 
Figure 10.1 illustrates the design of the entire CA-FE model. 


Temperature History within FE nodes 

< Linear Interpolation 

Temperature History in CA Model 

* i Materials Properties 

^ 

T 

Cellular Automaton Model 



Heterogeneous 

Dendrite Growth 

Grain Growth 

Nucleation 

Orientation 

Kinetics 


i 

Microstructure 



Figure 8.1 Flow chart of CA-FE model 
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Heterogeneous nucleation occurs nearly instantaneously at a characteristic undercooling. The 
locations and crystallographic orientation of the new nuclei are randomly chosen at the surface or 
in the liquid. The continuous nucleation distribution dn/dAT', which characterizes the 
relationship between undercooling and the grain density, is described by a Gaussian distribution 
both at the solid-liquid interfaceand in the bulk liquid (Eq. (8-1)). The parameters of these two 
distributions, including the maximum nucleation density ^ Ihe mean undercooling AT^, and 

the standard deviation of the grain density distribution AT^, can be obtained from experiments 
and grain size measurements. [24] 


^ ^ dn , 

Jq uAT Jq / 




AT^V^ 


exp 


1 [AT' - ATj^ 

2\ 


cLAT' 


( 8 - 1 ) 


Undercooling is the most important factor in the columnar and dendrite growth rate and grain 
size. The total undercooling of the dendritic tip consists of three parts, solute undercooling, 
thermal undercooling, and curvature undercooling. For most metallic alloys, the kinetic 
undercooling for atom attachment is small, so it is neglected [25]. The total undercooling can be 
calculated as follows: 


AT = mCo[l - A(Pc)] + 6^1 (Pt) + y 


( 8 - 2 ) 


where m is the liquidus slope; F is the Gibbs-Thomson coefficient; Cq is the solute concentration 
in the liquid far from the solid-liquid interface; and are the thermal and solutal Peclet 
numbers, respectively; k is the solute partition coefficient at the solid-liquid interface; A(Pc) 
equals [1 — (1 — 0^ is the unit thermal undercooling (= Ahfic); and R is the radius 

of the dendritic tip. 

For the laser deposition process, the rapid solidification condition corresponds to a high Peclet 
number at which the dendritic tip radius is given by (Eq. (8-3)): 

P nl/2 


R = 


[a*(mGc -G*)\ 


( 8 - 3 ) 


where a*, the marginal stability constant, approximately equals l/4n^ [26], and G* and G^ are 
the effective temperature gradient and concentration gradient, respectively. 
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8.2. Thermal History and Grain Morphology Model Result During Laser 


Remelting Process 


In this part of the experiment, no material was added on the substrate. The thermal history and 
grain morphology model in the molten pool domain during the laser remelting process appear in 
Figure 8.2 and Figure 8.3. Because the CA model requires a much finer grid than the FE 
temperature model, interpolation is necessary to obtain finer temperature information. The CA 
model domain corresponds to the black rectangle in Figure 8.2 
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Figure 8.2 Macroscopic FE simulated temperature field for laser remelting process on 
SS316. Temperatures are in K and shown at a) t=105ms, b) t=120ms, c) t=135ms, and d) 
t=150ms. The black rectangle indicates the section of thermal results used for the CA 
microscopic model. The rectangle is 200 pm wide and 70 pm high. 
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Figure 8.3 Microscopic CA simulation evolution of SS316 microstructure during 
solidification. Microstructure evolution is shown at a) t=105ms, b) t=120ms, c) t=135ms, and 
d) t=150ms. The ‘P and ‘J’ axes are the number of 2fim grid in the microscopic model, 
depicting the grain lengths in two dimensions. Color indicates individual grains. 
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Figure 8.3 indicates that equiaxed grains dominated the molten pool and shows the evolution of 
the solidification microstructure during the laser deposition process of SS316. The grid size for 
the macroscopic model was lOqm, and the macro time step was The grid size for the 

microscopic model was 2p.m, and the micro time step was 2 X As Figure 8.3 illustrates, 

the grains initially nucleated at the interface between the solid and liquid because the heat transfer 
at the interface occurred much faster than in the liquid. The different colors represent different 
crystallographic orientations of grains. Also, during grain growth, the grains always maintained 
their original orientations, which agreed with the physical mechanism. 

8.3. Solidification Grain Morphology Modeling Under Laser Deposition. 



Figure 8.4 Sample and cutting position scheme 
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(a) Front view 
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(b) Iso view 

Figure 8.5 Simulated temperature distribution during single-layer laser deposition process 


Figure 8.4 indicates the position of the deposition area on the substrate and the observation cross 
section. ABAQUS was used to simulate the temperature field of the substrate and deposited 
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material, including the deposition materials added on the substrate when the laser moved forward. 
A Fortran subroutine was written to simulate the moving laser’s heat source. The laser deposition 
of single layer Ti6A14V was conducted with the power of 750W, scanning speed of 600mm/min 
and powder delivery of 12g/min. The deposited material was added by activating elements in the 
deposition domain in a stepwise manner. The deposition temperature field and grain morphology 
were simulated first only in one layer. Figure 8.5 indicates the temperature field of the substrate 
and deposited material. 



T 

Figure 8.6 Grain morphology simulation of single-layer deposition of Ti6A14V. In the 
legend, ‘CLASS’ represents orientations of different grains. 

Figure 8.6 indicates that grain morphology simulation of single-layer deposition of Ti6A14V. The 
domain size in the CA model was 2mm x 0.4mm, and the simulation time was 100ms. The T’ 
and ‘J’ axes indicate the numbers of CA cells along the horizontal and vertical directions. The 
original temperature nodes obtained from ABAQUS were linearly interpolated because the 
element size in ABAQUS was 200|im. The size of each cell was lOpm in the CA model. 
Different colors indicate individual grains. The grain nucleation sites and growth rate were 
controlled by the temperature history obtained from the FEA model. For one-layer laser 
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deposition microstructure simulation, the grains were almost equiaxed. The grain size varied from 
30|im to 150|im. There are no columnar grains for one-layer laser deposition. 
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(b) Front view 

Figure 8.7 Thermal history for 25-layer Ti6A14V laser deposition 
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Figure 8.7 depicts the temperature field of the substrate and deposited material. ABAQUS was 
used to simulate the temperature field of the substrate and deposited material, including the 25- 
layer deposition materials added on the substrate when the laser moved forward and backward. 
The laser deposition of single layer Ti6A14V was conducted with the power of 750W, scanning 
speed of 200mm/min and powder delivery of 2g/min. 



8.8 Grain morphology modeling of 25-layer Ti6A14V laser deposition 

Figure 8.8 shows multiple layers of the Ti6A14V grain morphology under the laser deposition 
process. When more layers were deposited, columnar grains began to dominate, while equiaxed 
grains began disappearing and became the minority. The orientations of the columnar grains were 
almost perpendicular to the laser motion’s direction because the grains grew along thermal 
gradient direction. The domain size in the CA model was 2mm X 2mm. Based on this simulation 
result, the width of the columnar grains varied from lOOpm to 600|im. 
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8.4. Comparison with experimental results in the related literature 



Figure 8.9 Optical morphologies of DLF samples obtained: (a-c) at the bottom, center, and 

top locations 

Wu et al. [27] conducted a laser deposition experiment with a 1750-W C02 laser, showing the 
grain morphology results from using a range of laser powers, from 180 to 516W, and at a given 
scan speed of 300 mm/min, powder feed rate of 12 mm/min, and Z-increment of 0.3 mm. As 
shown in Figure 8.9, samples were fabricated at 264 W with a scan speed of 300 mm/min and 
powder feed rate of 12 g/min. The deposition height of each figure was 2.5mm. The total height 
of the deposition zone was 20mm. At the bottom, center, and top, the width of the columnar 
grains varied from 0.2 to 1 mm, 0.3 to 1.5 mm, and 0.2 to 4 mm, respectively. In the model, the 
columnar grain width varied from 100p.m to 600 p.m within a 2mm x 2mm domain. In the 
literature, the columnar grain width varies from 200p.m to 1mm within a 2.5mm x 
2.5mm domain. The current grain morphology modeling result is in the same order of magnitude 
as optical morphologies. 

8.5. Summary 

During the solidification process, the grain nucleation and growth are intensely controlled by the 
cooling rate of substrate and deposited zone. The thermal history can be simulated by the finite 
element method, which integrates with the cellular automaton method to predict the grain 
morphology. Grain morphology simulation of single-layer deposition of Ti6A14V indicates that 
equiaxed grains dominate in the deposited zone and no columnar grains can be found. Grain 
morphology modeling of 25-layer T16A14V laser deposition indicates that while more layers were 
deposited, columnar grains began to dominate and equiaxed grains began disappearing and 
became the minority. The orientations of the columnar grains were almost perpendicular to the 
laser motion’s direction. Compared to the similar studies in the literature, the multiple grain 
morphology modeling result is in the same order of magnitude as optical morphologies in the 
experiment. 
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In future work, the thermodynamic properties and parameters of compositions and phases, which 
can be obtained from a commercial database such as JmatPro, can be included to control grain 
nucleation and growth. This will improve the accuracy of the grain morphology simulation result. 
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9. Microstructure Observation and Simulation of Ti6A14V Single-Wall 

Deposition 


9.1. Introduction 

The nature of a/p transformation in Ti alloys and the complex thermal history in the laser 
deposition process causes some special morphology and microstructure features to appear in the 
deposits. 

Different Ti alloy microstructures are characterized by the size, shape, and distribution of two 
phases, a (hep, below P-transus) and P (bcc). P transforms into various products in the a/p region. 
The microstructures appearing in LMD include colony a, plate-like a, grain boundary a, 
Widmanstatten a, and martensite a’. The latter is produced upon rapid cooling (>16.7 K/s [28], 
>410 K/s [29]) above the martensite start temperature and can transform into a very fine lamellar 
a/p microstructure between 700-850 °C. An intermediate cooling rate produces Widmanstatten a, 
while a comparatively slow cooling rate results in colony a. The labels Widmanstatten a, 
basketweave a, and acicular are often used to refer to the same structure. 



Figure 9.1 Various typical microstructures appearing in Ti alloy laser deposits [30]. (a) 
shows colony, grain boundary and basketweave a; (b) shows basketweave a, martensite a’ 


and prior p-grain boundary 

An obvious transformation trend exists from colony a and basketweave a to martensite a’ as the 
cooling rate increases at P-transus. However, the critical transition rate is not clear. The 
relationship between the a-lath thickness, the size of colony a, the microstructure percentage, and 
the cooling rate and other temperature information is the key to predicting the microstructure in 
laser deposits. 
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9.2. Experiment 

To study the evolution of the microstructure during the laser metal deposition process, a 
deposition experiment was conducted using a diode laser (Nuvonyx ISL-IOOOM, 808nm) with a 
maximum power of 1 kW. The laser radiation energy density was evenly distributed on a circle 
with a 2 mm diameter. The laser beam head, powder feeder nozzle, and substrate were fixed on a 
CNC machine to create the deposition path. 

Thin walls were built according to the parameters listed in Table 9.1. The powder and substrate 
were both T16A14V. The height of the deposit was 2 and 4.8 mm. 

Table 9.1: Experimental parameters 


Number of 
layers 

Layer thickness 
(^m) 

Laser power 
(W) 

Scan speed 
(mm/min) 

Scan length 
(mm) 

25,60 

70-80 

750 

200 

10 


The deposits were cut through the middle cross section perpendicular to the direction in which the 
laser moved, and through the central cross section parallel with laser’s direction, by wire 
electrical discharge machining (wire EDM), as shown in Figure 9.2. The specimens were 
mounted, ground, polished, and etched according to standard metallography methods. An optical 
microscope was used to observe the morphology. 




y..L,K 

^vLx 

(a) First cross section perpendicular to the (b) Second cross section parallel with 
laser’s direction laser’s direction 

Figure 9.2: Cross section position for optical microscope microstructure observation, (a) 
First cross section perpendicular to the laser moving direction; (b) Second cross section 

parallel with laser moving direction. 
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9.3. Simulation 


The finite element analysis (FEA) thermal simulation was accomplished using Abaqus 6.13. The 
model parameters were the same as in the experimental setup, as shown in Table 9.1. The laser 
radiation was modeled as a moving heat flux condition via a user subroutine, and new material 
was added into the deposit via an activation/deactivation mechanism in Abaqus. 



Figure 9.3: Simulated temperature distribution during 25 layers Ti6A14V alloy deposition 

at the first layer. 

Because of computer powder and time limitations, only 25 layers of the sample were simulated. 
Figure 9.3 shows the temperature distribution result during the deposition of the first layer. 


9.4. Result and Analysis 


The microstructure observations of 25 layers at several distances from the top surface along the 
central vertical line of the deposit appear in Figure 9.4. Clearly, all of the microstructures were 
typical a’ martensite. Fong, orthogonally oriented martensitic plates are clear in these regions. At 
the bottom, some indication of basketweave a appears on the right side of the image. However, 
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Figure 9.4: Microstructures at different locations of the central vertical line from the top to 
the bottom of the 25-layer deposit, (a) top; (b) 0.25 mm below top; (c) 0.75 mm below top; 
(d) 1.25 mm below top; (e) 1.75 mm below top; (f) 2 mm below top (bottom) 


The microstructure results of the 60-layer sample appear in Figure 9.5. At the top of the middle 
cross section, the a’ martensite formed inside the prior p grain. At 3 mm from the surface, 
orthogonal acicular morphology still existed. However, compared with the martensite in the top 
area, the acicular lath was thicker, and long plates broke into smaller ones. The morphology 
changed at 4.2 mm from the surface (Figure 9.5c). Two areas having different microstructures 
can be discerned. Area A was still a’ martensite, while Area B consisted of basketweave a laths, 
and no a’ orthogonal grid feature appeared. At the bottom of the deposit, the morphology was 
typical basketweave a. 
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Figure 9.5: Microstructures at different locations of the central vertical line from the top 
surface of the 60-layer deposit, (a) top; (b) 3 mm below top; (c) 4.2 mm below top; (d) area 
A marked in (c); (e) area B marked in (c); (f) 4.8 mm below top (bottom) 


The simulated temperature history of the top and bottom locations in the 25-layer deposit appears 
in Figure 9.6. The cooling rates at the P transus temperature of Ti6A14V (1243 K [31]) are also 
indicated in the figure. The temperature of the top region during the last scan approximately 
reached the boiling point. The bottom region also heated to above the P transus during every laser 
layer scan. Therefore, the a/p microstructure at the bottom region formed during the last laser 
scan. 


65 







Figure 9.6: Temperature history at top and bottom locations of 25 layers simulated 


deposition result. Cooling rate p-transus point (1243K) is indicated. 

The cooling rate from the p phase region to the a region determines the formed a morphology 
[32]. Some studies have concluded that cooling rates greater than 18 or 20 K/s lead to the 
development of a’ martensite [3 3] [34]. For the 25-layer sample, the cooling rate distribution 
along the central vertical line appears in Figure 9.7. Every position in the center may experiences 
several reheating and cooling processes and the cooling rate mentioned here refers to the last time 
that determines the formed microstructure. Cooling rate increases from the bottom to top region. 
This is due to although the cooling rate at bottom region at the first time maybe very large, when 
cooling from Tp at the last time temperature is much more evenly distributed due to thermal 
diffusion. The cooling rate results and 20 K/s criteria verify that a’ martensite dominated the 
microstructure of the 25-layer sample, as depicted in Figure 9.4. 
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Figure 9.7: Cooling rate at 1243 K distribution along the central vertical line of 25 layers 

simulated deposition result. 


9.5. Conclusion 


The nature of the a/p transformation in Ti alloys and the complex thermal history from the laser 
deposition process cause some special morphology and microstructure features to appear in the 
deposits. The microstructures of 25- and 60-layer samples exhibited different distributions. The 
former featured a’ martensite, whereas the latter transitioned from a’ martensite to basketweave 
a. While the cooling rate of the 25-layer deposit at 1243 K conformed to results published in the 
literature, the transition cooling rate could be determined only if the temperature simulation of the 
60-layer sample was available. 

Further experiments and simulations on larger deposits need to be conducted to reveal different 
microstructures of Ti6A14V in the laser metal deposition process. 
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10. Incorporation of Infrared Camera to Monitor Laser Metal Deposition 

Process 


10.1. Introduction 

A FLIR A615 camera, an industrial automation infrared (IR) camera with sensitivity in the range 
of 7.5-13 microns, was used to perform the following studies. The capabilities of the camera are 
listed in Table 10.1. 

Table 10,1 Specifications of the IR camera 


Feature 

Specification 

Spatial resolution 

0.69 mrad 

Focal length 

25.4 mm 

F-number 

1 

Imaging frequency 

12.5 Hz to 200 Hz 

Image resolution 

640x480, windowing at high freq. 

Temperature measured 

3 ranges, -50 C to 2000 C ( e=l ) 

Detector 

Uncooled bolometer 

Detector time constant 

8 ms (typical) 

Spectrum sensitive to 

7.5 to 13 microns 


68 


Initial efforts focused on realizing the setup requirements and possible calibration steps to use the 
IR camera as a process monitoring tool. Later, the camera was used to see the temperature trends 
during deposition and the size of the mushy zone while using the LAMP lab’s feedback. System 
and acquired data were processed to identify the melt pool or the mushy zone and solidus regions 
in the deposit. 

10.2. Setup and Calibration 

An IR camera measures the radiative power of a hot body. This power is later converted into a 
temperature using lookup tables; hence, meeting the established parameters is crucial for 
acquiring accurate temperature readings. The laser metal deposition process is very dynamic, and 
several possibilities can result in significant and irregular changes in the input parameters 
required for the temperature conversion. The most pertinent of these parameters is emissivity; for 
accurate measurement, the spectral and thermal dependencies of the emissivity must be identified 
and accounted for. For the initial set of experiments, the thermal and spectral dependencies were 
assumed to be negligible. 

The camera performs only gray body measurements, meaning the acquired data can be processed 
with a single emissivity (unless a pixel-by-pixel emissivity map is provided). The associated 
literature states that the emissivity of metal oxide is higher than the respective metal, and 
emissive properties are favored by higher temperatures. Therefore, to rationalize the assumption 
of emissivity’ s temperature independence, the acquisition was performed in an open atmosphere 
that resulted in significant surface oxidation (considering the oxide’s emissivity at room 
temperature was approximately 0.9). 

10.3. Experimentation and Results 

The deposition of stainless steel 316 was monitored using the IR camera. The energy 
management system was tested and monitored to determine its effect. The emissivity was set to 
0.95, as arrived at by averaging the emissivity values of the oxides of constituent elements. A 
capture rate of 200 fps was set for thermal history acquisition, and other parameters, such as the 
atmospheric temperature, humidity, atmospheric reflective temperature, and atmospheric 
transmission, were set as measured. 

10.3.1. Energy Management System 

The energy management system (EMS) manages the amount of energy according to a certain set 
signature, i.e., it is set to monitor the radiation and maintain the energy below a set point by 
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modulating the laser power. The experiment was performed for three different set points, and the 
power logs in Figure 10.1 and Figure 10.2 illustrate the effect of each set point. 


Po^r V&riaiiork with Time 
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Power Vari alien wiih Time 
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Figure 10.1 Power logs for set points 1 and 2 
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Figure 10.2 Power log for set point 3 

The cutoff increased in value from set point 1 through 3, i.e., the allowance for the amount of 
energy contained in the deposit at a particular instance increased. The power modulation tended 
to undergo an exponential decay of sorts as the deposition progressed. 

10.3.2. High-Temperature Region 

The minimum cutoff temperature for this region was set as the solidus temperature of SS 304. 
The assumptions imposed during parameter selection create an unknown amount of error in the 
temperature readings acquired. Nonetheless, this experiment was instrumental in identifying the 
size trends of this high-temperature region during deposition. 

10.3.3. Deposition Visualization 

Using the iron color palette, a false color image was used to visualize the temperature data. The 
region above the minimum cutoff temperature also was visualized through the deposition. 



Figure 10.3 False color images used to visualize the deposition and highlight the high- 
temperature region at an instance during deposition 
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For each of the noted set points, the average powder feed rate was changed, and the size trends of 
the high-temperature region were plotted as the deposition progressed. Figure 10.4 illustrates 
these trends. 
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Figure 10.4 Size trends of the high-temperature region through deposition for variable 

powder feed rates and set points 
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The optimum solution for ensuring homogeneity throughout the deposition process would be one 
that saturates the size of the high-temperature region. The current experiment, through trial and 
error, was able to successfully realize such saturation, thereby identifying the possible optimum 
setup for homogenous deposition. Though not a smooth saturation, the size of the high- 
temperature region stabilized around a certain range for a set of parameters. 

10.4. Image Processing 

The associated literature [35] states that a drop in emissivity occurs when a metal undergoes 
phase transformation, i.e., when solid is melted into liquid, the emissivity decreases. As the IR 
camera reads the radiative power, it senses a decrease in the temperature when such a 
transformation occurs (as the temperature is analyzed using a single emissivity value). The band 
radiance pertinent to the camera is listed in Table 10.2. 

Table 10.2 Band radiance pertinent to the camera 


SS316 

Temperature 

Emissivity[36] 

Band radiance 
(7.5 to 13 micron) 

Liquidus, 1400 C 

0.3 

4745 W/sq.m/sr 

Solidus, 1377 C (oxidized) 

0.9 

1582 W/sq.m/sr 


This proves that the melt pool appears colder than it actually is. Exploiting this phenomenon, 
edge detection techniques were chosen to identify the melt pool or the mushy zone in the 
thermographs. The change in emissivity is illustrated in Figure 10.5. 



^ solid ^ ^ melt pool/mushy zone ^ ^ ambient zone 


Figure 10.5 Boundaries and emissivity change trend during deposition 

1 0.4. 1 . Laplace Edge Detection 

The Laplace edge detection technique was chosen to identify the transitions in the deposit. 
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Figure 10.6 Steps in Laplace edge detection [37] 

The Laplace edge detection technique targets the change in intensity of the signal and identifies 
the point of transition using the zero crossing points from the second-degree derivative of the 
signal. The points of interest or the edges are identified using a threshold variance as the decision 
criterion in step 3 (from Figure 10.6) [38]. 

10.4.2. Methodology 

Steps involved in image processing: 

1) Moving median 

2) Gaussian blur 

3) Laplacian transform 

4) Find zero crossings 

5) Clarify image 

6) Detect edges 

7) Identify and process region around melt pool 

8) Mark melt pool and solidus region 

These steps were followed to remove noise, smooth, identify edges, clarify the image, and mark 
the region of interest. 

The outcomes of implementing these steps appear in Figure 10.7. 

□ IDQ 200 30Q AOn SDO 60D 

(a) Moving median 
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Figure 10.7 Initial image processing steps 

The edges found, as shown in Figure 10.7, do not include the melt pool because the variation 
from solid to liquid was insignificant compared to the variation from the deposit to ambience. In 
order to identify the melt pool, a more local search had to be performed within the deposit. After 
further image clarification and edge detection, the final outcome was arrived at, as depicted in 
Figure 10.8. 
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Figure 10.8 Deposit (sky blue), solidus region (yellow), and mushy zone (red) 

10.5. Conclusions 

An infrared camera was incorporated as a process monitoring tool to identify the solidus and 
mushy zones during deposition. The images were successfully processed to identify these regions. 

10.6. Future work 

The relationship between the temperatures in the regions of interest needs to be established, and 
further emphasis must be placed on decoding the constituents of the mushy zone. These 
explorative efforts need to be employed during deposition to identify various phenomena during 
solidification. 
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phenomena associated occur in multiphysics and multiscale. 
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